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ABSTRACT 


The theoretical response of the earth to elastic 
waves is studied by evaluating an analytic solution of 
the elastodynamic equations of motion. 

Spherical earth models with radial ee Of 
the elastic parameters are studied by applying an earth 
flattening transformation to the model, and approximating 
the flattened model by a stack of plane homogeneous 
layers. 

the exact Cagniard-de Hoop theory is’ used to find 
the response for each important ray in the ray expansion 
of the complete response in the plane homogeneous layers. 

Numerical long period seismograms are presented 
for SH, P, and SV waves near the core geometric shadow 
edge, and deep into the shadow, in a homogeneous earth 
model. Comparison of the seismograms and their spectral 
amplitudes with frequency domain wave solutions of the 
same problem confirms the applicability of this time 
domain ray method to modelling studies of the lower man- 
tle and the core-mantle interface. Previous ray methods 
are not applicable near the core shadow edge. 

Numerical seismograms at ranges of 84° to 96° 
have been computed for a model by Birch (1964), and for 
other structures derived from the Birch model by 


Kanasewich et al (1973) to investigate a region of 
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lateral heterogeneity in the lower mantle discovered by 
seismic array observations. The Cagniard-de Hoop ray 
method is the most appropriate method for constructing 
numerical seismograms near the core shadow edge for 
models with thin resonating layers. A comparison of 
the results with observational data shows good quali- 
titative agreement, and suggests that the Birch model 
needs alterations near the core-mantle interface. 

There is a need for a two and three dimensional modell- 
ing Capability to fully “resolve:fine ‘structure in the 


lower mantle. 
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CHAPTER. 1 


INTRODUCTION 


Lok Numerical Seismograms 


The greater part of our knowledge of the structure 
of the interior of the earth has come from seismology, 
the study of the transmission of elastic waves from 
earthquakes and explosions. The major discontinuities 
within the earth were identified, and the average depth 
dependence of the elastic parameters grossly determined 
prior to 1940, mainly through the use of travel times of 
identifiable types of seismic pulses. Mohorovicic, in 
1909, identified the crust-mantle interface. Gutenberg 
in 1914 determined the depth to the central fluid core 
postulated in 1906 by Oldham, and Miss Lehmann in 1936 
obtained evidence of the existence of a solid inner core. 

The Jeffreys-Bullen travel time curves of 1940) 
give arrival times for all geometric types of elastic 
wave to accuracies of a few seconds and are widely used 
by Garth scientists throughout the world to identify 
seismic arrivals. 

With the introduction of more reliable instrumen- 
tation and improved spatial distribution of observations 
by the Standardized Stations of the USCGS worldwide net- 


work and by the LRSM Stations, it is now feasible to 
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determine earth structure in much finer detail by 
comparing elastic wave amplitude data with amplitudes 
predicted by solving equations of motion numerically 
on the large digital computers now available. Numeri- 
cal seismograms produced in this way are most useful 
in the interpretation of data where simple geometric 
tay theory (Bullen, 1963, p.109) as inappropriate, ive. 
at any discontinuity in the travel time curves. Dif- 
fracted or non-geometric signals are observed in the 
geometric shadow of the earth's liquid outer core and 
solid inner core (Bolt, 1962), in the shadow of the 
upper mantle low velocity zone, and at the end points 
Gf any ‘triplications due to high velocity gradients in 
the mantle (Johnson, 1967). 

Numerical seismograms are computed by choosing an 
appropriate earth model defined by its elastic parameters, 
introducing a source of elastic waves, and solving the 
elastodynamic equations of motion to obtain the response 
of the earth model to the elastic waves at any desired 
place and time. 

The vector, second-order, partial differential 
equation of motion is most easily solved through a 
Laplace or Fourier transform method. This leads 
naturally to two types of solution. Wave solutions 
solve the time transformed equations for the total 


spectral response at a particular frequency. A band- 
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limited frequency spectrum can be constructed by solving 
at various frequencies. The time domain response at all 
times is obtainable by an inverse Fourier transform. 
Normal mode solutions (Alterman and Aboudi, 1969), 
matrix solutions (Haskell, 1953; Randall, 1967), and 

the numerical technique of Chapman and Phinney (1970) 
are in this category. 

Alternatively, all inverse transformations are 
often included in the solution process. This produces 
directly the time domain solution containing all fre- 
quencies in a time-limited window. The asymptotic ray theory 
(Alekseev et al,1961)and the exact generalised ray theory 
(Cagniard, 1939; Pekeris et al, 1965; Helmberger, 1968) 
are examples of time domain solutions. Time domain 
solutions using Cagniard-de Hoop exact ray theory are 
discussed in Chapters 2 and 3. 

The structure of the core-mantle interface and 
of stherBud len region D™) (Bullen, 1963, p#@222) above 
thei interface is still uncertain. Chapter 6 discusses 
a study illustrating the feasibility of examining the 
structure of this region by the use of Cagniard-de Hoop 
seismograms for the directed and reflected waves in the 
illuminated region and the diffracted signal in the 
geometric shadow of the core. The results are compared 


with published wave solutions for the same problem. 
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Chapter 7 presents a study modelling heterogeneities 
in the lower mantle under Hawaii suggested by the 
seismometer array studies at the University of Alberta 


(Kanasewich et al, 1972). 


1.2 Wave Solutions 


Solutions of the equations of motion transformed 
over time and horizontal distance in an inhomogeneous, 
radially symmetric earth model can be expressed as a 
sum of spherical harmonics at each frequency to obtain 
normal mode wave solutions. The solution is thus ex- 
pressed in terms of the natural free oscillations of 
the earth. This method works well for surface waves, 
but many terms are necessary to include frequencies 
high enough to resolve P and SV body waves. 

The transformed solution of the equations of 
motion in a homogeneous earth with a fluid core can also 
be written as a sum of spherical Bessel functions at 
each frequency to obtain a rainbow or Debye ray expan- 
sion (Debye, 1908) equivalent to a partial wave expansion 
in scattering theory. 

This method becomes difficult when the earth 
model is inhomogeneous. An inhomogeneous earth model 
can be approximated by a series of concentric homogen- 


eous shells. Within each layer the response can be 
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obtained in “terms “of “spherical Hankel “functions seat 
each interface, the normal displacement and traction 
must be continuous. These boundary conditions create 

a matrix equation to be solved for each frequency 
component. Haskell (1953) used this method for surface 
wave studies, and Chapman (1969) reformulated the tech- 
nique to study body shear waves. 

An alternative form of wave solution for inhomo- 
geneous models was introduced by Chapman and Phinney 
(1970). Using generalised displacement potentials 
(Richards, 1971), an approximate solution of the trans- 
formed equations of motion can be obtained using WKBJ 
asymptotic solutions (Chapman and Phinney, 1972). In 
the region of the: turning point of the wave, where the 
WKBJ approximation is invalid, the transformed equations 
of motion are solved numerically instead of analytically. 
Each spectral component of the direct wave and of the 
core reflection in the illuminated region is found by 
a contour integration over two saddle points in a com- 
plex wavenumber space. Spectral components of the 
diffracted signal in the geometric core shadow are 
found by evaluating the residues at a series of poles 
in the same space (Chapman and Phinney, 1972). 

This process is carried out for typically sixteen 


equally-spaced frequencies. Up to fifteen intermediate 
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frequencies can be reasonably interpolated between each 
pair of calculated components, to obtain a band-limited 
frequency spectrum. An inverse Fourier transform pro- 


duces the time domain signal. 


alae) Ray Solutions 


Time domain solutions of the equations of motion 
can be expressed as a sum of rays. A ray, represented 
by a continuous line from the source to the receiver, 
is one possible path for energy travelling to the 
receiver. In a non-dispersive medium, the ray path 
is at all points orthogonal to the wavefronts of the 
disturbance emitted by the source. 

Lamb (1904) first studied the time domain resporse 
of a homogeneous half-space to an imbedded disturbance, 
and Cagniard (1939) obtained a time domain solution for 
the single ray by applying Laplace and Fourier transforms 
to the equations of motion to obtain the transformed 
solution, and analytically solving the inverse transform 
integrals. 

Pekeris et al (1965) extended the method to the 
problem of a solid layer on a solid half-space and were 
able to approximate the complete time domain solution 
over a wide range of distances by summing the response 


of a few simple rays. 
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The exact solution for an arbitrary ray in a 
stack of homogeneous layers was discussed by several 
authors (Podyapolski, 1959; Helmberger, 1968). 

Numerical seismograms using this exact ray 
theory in multi-layered models were still impractical 
due to the inverse transforms. De Hoop (1960) intro- 
duced a simplification to the solution of the half-space 
problem (Lamb's problem) so that two of the inverse 
transform integrals were eliminated. Gilbert (1963) and 
Helmberger (1968) applied de Hoop's method to the multi- 
layered model, and realistic numerical seismograms became 
practical on digital computers. Helmberger (1968) inves- 
tigated crustal layering and the crust-mantie transition 
in the Bering Sea using this Cagniard-de Hoop method. 
Helmberger and Wiggins (1971) investigated upper mantle 
structure in the United States by synthesizing seismo- 
grams for NTS nuclear explosions. 

It then became feasible to approximate models with 
continuous gradients of parameters by means cof stacks of 
thin homogeneous layers. Since no individual ray reflec- 
ted by a discontinuity in the layered model corresponds 
to the wave refracted by the velocity gradient in the 
original model, it was necessary to develop a ray 
expansion whose sum in the layered model could be shown 


to be equivalent to the wave in the original continuous 
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model. Hron (1971) has classified all possible rays 

in a model of plane homogeneous layers. An expression 
for a complete ray expansion can be derived from this 
(Hron and Chapman, personal communication). Cisternas 
et al (1973, in press) have shown that the ray expansion 
solution in the layered model is the complete solution 
in that model. Muller (1970) has outlined a procedure 
to determine the sufficient nUnes Sir layers valde as sut— 
ficiently complete ray expansion to ensure equivalence 
of the layered model response and the response of the 
original continuous model. It is now possible to cal- 
culate approximate numerical seismograms for any one 
dimensional structural model by using a truncated series 
of generalised rays. 

Exact ray theory is not directly applicable to 
teleseismic studies because of earth curvature. Gilbert 
and Helmberger (1972) used an approximate form of the 
exact ray theory to obtain short period pulses reflected 
in a sphere of homogeneous layers. The method is not 
applicable to either direct waves in the mantle or the 
diffracted signal in the geometric shadow of the core. 

A more useful approach is to apply an appropriate 
earth-flattening transformation to the elastic parameters 
of the spherical earth model to construct a flat model 


whose response to an imbedded disturbance closely 
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approximates that of the original spherical earth model 
at all times and distances of interest. The exact ray 
theory can then be used to study the flat model. 

Muller (1970) has computed numerical seismograms for 

an upper mantle model of Julian and Anderson (1968) 
after using such an earth-flattening transformation. 

Earth-flattening transformations have previously 
been used in surface wave studies. Biswas and Knopoff 
(1970) presented a transformation that is exact for Love 
waves. Biswas (1972) and Schwab and Knopoff (1972) used 
another form of power law transformation of elastic para- 
meters to study Rayleigh waves. Helmberger (1972) has 
used a power law earth-flattening transformation for 
body wave studies. 

Chapman (1973, in press) has determined the opti- 
mum form of power law transformation for various condi- 
tions. In a solid elastic model, the Biswas and Knopoff 
(1970) transformation is optimum for all body waves. 

It has been used in this study. 

Another form of ray theory applicable to mere 
general models than the exact ray neces has also been 
discussed (Alekseev et al, 1961). The high frequency 
response for each ray in the ray summation is approxi- 
mated by the leading terms in an asymptotic series. 
Higher terms.are needed for rays reflected “at either 


near-vertical or near-grazing incidence (Podyapolski, 
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1966). The method is not applicable to rays near the 
critical range where the reflected wave and the head 
wave arrive nearly simultaneously. At all other dis- 
tances, this asymptotic ray theory response can be 
obtained much more easily than the exact ray theory 
response. Asymptotic ray theory is well suited to the 
study of geological structure, and has been used to 
interpret seismic reflection data in Alberta (Hron and 
Kanasewich, 1972). 

Although it was not used in this study, asympto- 
tic ray theory could be combined with exact ray theory 
for critical and grazing rays, to study more general 
earth models efficiently. 

This study demonstrates the feasibility of using 
the Cagniard-de Hoop method and exact ray theory with 
the Biswas and Knopoff earth-flattening transformation 
to study diffracted waves in the geometric shadow of 
thesearnth\isiccoreby duplicating mesults) of diffraction 
studies by Chapman and Phinney (1972) who used a wave 


method. 


1.4 The Core-Mantle Interface 


Detailed knowledge of the nature of the core- 
Mantle “interface is necéssary for a complete understand- 
ing of the formation of the earth and the generation of 


its magnetic field. 
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The history of the formation of the liquid core 
may be resolved by the determination of the thickness 
of the core-mantle transition, atid the nature of the 
chemical and phase changes, and the temperature. 
gradient across it. Current evidence suggests the 
exrstence ofa thin, well-detined intertace, overlain 
bya transitivon Layer 30 to 150 km. thick, of lower 
seismic velocity, higher density and higher temperature 
gradient than the mantle above (Cox and Cain, 1972). 

Small variations in the length of the day, of 
che OLrderaor Kiger seconds seasonally and fo seconds 
over periods of years have been measured accurately 
SinecesthewinerOGuction OL atomic clocks. The vexplana— 
tion of these variations appears to require a transfer 
of angular momentum across the core-mantle interface. 
Hide and Malin (1971) suggested topographic bumps of 
the order of 10 km. over several hundred kilometres 
Might provide a mechanical coupling mechanism. 

The search for these bumps indicates that undu- 
lations Larger than-2 km. do not exist (Bolt, 1973), 
due to consistent arrival times of waves reflected up 
to seven times inside the core. The search has, however, 
produced evidence of important lateral heterogeneities 
in the mantle and core that may provide the necessary 


coupling by electromagnetic processes. 
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Studies of the non-dipole geomagnetic field and 
reversals in geomagnetic polarity also indicate inho- 
mogeneities in electromagnetic properties and irregu- 
larities in convection in the fluid core near the core- 
mantle interface. More detailed knowledge of this 
region of the earth is needed to understand the dynamics 
of the geomagnetic field. 

The hypothesis of "hot spots" with a scale length 
of a few hundred kilometres at the base of the mantle 
has come from three different types of investigation. 
Cox and Doell (1964) suggested that the irregularities 
in the geomagnetic field are controlled by lateral tem- 
perature variations at the core-mantle boundary. 

Morgan (1972) extended Wilson's (1963) ideas on island 
chain formation in the interior ,of a lithospheric plate 
to suggest that the lava of island chains such as Hawaii 
originates near the core-mantle boundary in localized 
hot areas, and rises in the form of a narrow plume in 
the near-stagnant centre of a convection cell. 

Seismic array studies (Kanasewich et al, 1972; 
Davies and Sheppard, 1972) have revealed evidence of 
lateral velocity variations at the base of the mantle 
under Hawaii. These anomalous velocities may be indi- 
cative of temperature variations. Numerical seismo- 
grams for a velocity model of the Hawaii area 


(Kanasewich et al, 1973) have been computed in 


Chapter 7. 
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Diffraction of elastic wave energy around the 
earth's fluid core was recognized by Wiechart within 
a few years of the discovery of the core. Jeffreys in 
1938 ‘discussed observations of diffracted waves at 
distances of up to 150° from the sotirce earthquake, 

i.e. more than 50° into the geometric shadow. 

Because this diffracted energy travels a large 
distance near the core-mantle boundary, the diffracted 
signal amplitudes are sensitive to the physical proper- 
ties of both the mantle and the core at their interface. 

Diffraction of waves by a sphere has been studied 
by Van der Pol and Bremmer (1937) for radio waves, and 
by scholte (1956), and Duwalo and Jacobs (1959) for 
high frequency elastic waves in a solid mantle surround- 
ing a fluid core. Nussenzveig (1965) gave asymptotic 
Solutions tor all regions £or acoustic ditiraction by an 
impenetrable sphere. 

Phinney and Alexander (1966) used a matrix wave 
method to examine the effects on the spectral amplitudes 
of diffracted signals due to layering and a low velocity 
zone at the base of the mantle. Alexander and Phinney 
(1966) have then interpreted spectral ratios of dif- 
fracted P waves for several large earthquakes and found 
evidence for lateral variations of lower mantle proper- 


ties, and for a low velocity: layer. 


as vevey inet ty cae bevaucabe Seer 

Alailphsuis stakisn Sit Rc OMORT of Gu 50 ANNES = 

waters Sixtecosp “sds-odni "02 medg és0m .osk 

seael © sievers yoians Sascsa2 tid aids saweosd ; 
Bedositirh sad yvrebruvced! Siinga-s2es sie iad vonasedb 
~fegomq Issieytg sft of svistense sis eablatiqns faaple 
 ~bositesni sist 48 ot05o Sds\ bas slinsn sid Micd Yo ests 


y te 


bethess ased sat Ssi1erge 4 yd eayew to NOitosexI3iG 
bee jaevew oibes sod ((COL) aamesrvd Soe icd =6b as yd 
wO2 (201) edooet bos clawed Sas . (deel) stfodse Yd - 
-Bauossve Sitdem Difoe 6 ot esvew piszsle yonoupss? Apid 
sisodaqmyas Svse (2081) pisv¥sassay% .aton Sinl? & pak 
ie Yd aoljosittib altstoos +02 stoipox (ls +0? sadbisuioe 
-atadge aideutanegal 
SvAW xEnIem & beey, (coli) ishaexeiA bas ysanidd 
ashusiiqus i6idosqe: orld Pere) efielis oi) siimexs os bodgem 
dara Alli hath Seu ini had een ieee 
eck 314 Loltnem ony 20 shed Sit 26 snes | 
ss Seeioninrmate: nt? sved (838) 
7 isegesialt -esvaw 2 bedoga?! 
be: veces Ges is spaahive 


14 


Phinney and Cathles (1969) have computed spectral 
amplitudes for P waves diffracted by a homogeneous core 
in a homogeneous mantle. Chapman and Phinney (1970) 
and Richards (1970) extended the method to other types 
CEmWave srs 

Knopoff and Gilbert (1961) calculated the asymp- 
totic diffracted pulse deep in the shadow by analytically 
performing inverse transforms for the residue of the 
first drriraction pole déscribed: imeSection ri .i2). Their 
results compare favourably with the observed diffracted 
wave amplitudes published by Gutenberg (1960). 

The pulse reflected from the core in the illumi- 
nated region can be computed by the approximate ray 
method of Gilbert and Helmberger (1972), but ray method 
solutions have not been successfully applied to signals 


near the geometric shadow previous to this study. 
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CHAPTER. 2 


GENERAL THEORY 


Zeb Squat lOns Of Motion 


The purpose of synthetic seismograms is to present 
the expected displacement u(t,r) at a point r as a func- 
tion of time t, resulting from a disturbance at another 
point Xo in an elastic medium. 


The equation of motion at the point r is 


a7u(t,r) 
Sore ree oe mS mera (2.1) 
at ss 


Wi Che py eoeanderethesmatertaledensity, the stress tensor 
and the body force per unit volume. 

When the medium is isotropic and homogeneous the 
Stress LeCnsOuy can abeswrhittenean ‘terms of u(t7n), where 


A and uw are the Lamé elastic parameters 
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In cylindrical coordinates (r,6,z), the components E45 


of the strain tensor are 
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Abtere thes substitutions (2.2) and. (2.3), Ene equation 


of motion (2.1) becomes 


2 
dg u(t,r) 5 
OP tN et AV UE, Bt) a i ee ey 
at 


The displacement and body force can be expressed 


in terms of potentials (Hook, 1961) 


u(t,r) = Vo + VxVxy, + VxP. 

Pz (t,¥)= (0, O, Y (t,r)) 

py (t,x) = (0, 0, Wy, (t,x)) (2-15) 
E£(tjxr) = (A+ 2n)V L + plVxvxM, + VxM] 

M(t,r) = (0, 0, M,(t,x)) 


M,, (t,x) =(0, 0, M,,(t,x)) ° 


The equation of motion (2.4) separates into three scalar 


equations for the potentials (Hook, 1961) 
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Equation (2.6) describes compressional waves (P waves) 


travelling at the compressional velocity a. 


Q = det 2 (2.9) 


The horizontal and vertically polarized shear 
waves are described by (2.7) and (2.8). The shear 


VeLOCZty USsyer 
B = Vu/p (210) 


The complete displacement response in an isotropic 
homogeneous elastic half space is obtained by solving 
C206 ) F257)" andat2k 8)” fortthes potentialss $4 Vay and Vy 
applying boundary conditions vat. theofree) surface, , then 
Seiving (295)' for u(t)r) .eHowever a@most mediatofeanterest 
are inhomogeneous. Richards (1971) has shown how the 
equation of motion (2.1) can be separated using potentials 
that represent P, SV, and SH waves in a radially inhomo- 
geneous sphere. However, the method is not applicable to 
the Cagniard-de Hoop formulation applied in this study. 


This study used plane, parallel homogeneous elastic layers 


(v.48) 
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with changes of density and seismic velocity at the inter- 
faces. The complete solution of the equation of motion 
(2.1) for such a structure may be obtained by the use of 
Fourier transforms and Haskell matrices (Haskell, 1953), 
satisfying boundary conditions on stress and displacement 
at each interface (see Alexander and Phinney, 1966 for the 
application of the method to body waves). The inverse 
transforms to synthesize a realistic causal pulse are 
laborious (Chapman and Phinney, 1972). Another matrix 
wave method by Knopoff (1964) and modified by Randall 
(1967), is also efficient to obtain the complete response 
to plane wave sources, but also requires a summation of 


plane wave frequency solutions for a realistic source. 


Zee) “The Ray Expansion 


Instead of summing the exact plane wave response 
at each frequency for all frequencies of interest, it 
is practical to calculate the response for energy arriv- 
ing at a given moment regardless of frequency, then sum 
over the time window of interest. This is a time domain 
method instead of a wave method. To achieve this, the 
Laplace transform of the response can be written as an 
infinite expansion of generalised rays in (s, p, &, 2) 
space, where s gives frequency, p and & are horizontal 
wavenumbers and z is the vertical coordinate. The 
response due to each ray is readily obtained and the 


total response is then found by summing rays. The 
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development of the ray expansion was not a part of this 
study. It has often been assumed in published ray stu- 
dies (Helmberger 1968, Muller 1970) and has been proved 
to be a complete solution by Cisternas et al (in press). 

A ray can be shown diagramatically in coordinate 
Space by a single connected series of line segments from 
the source™to™ the receivers "The ray patil! is straight 
within each layer, and represents the propagation direction 
of a wave through the layer. The ray path is at every 
point the normal to the wavefront for all plane wave com- 
ponents in the source spectrum. At each interface encoun- 
tered, the ray is either reflected or transmitted at the 
appropriate angle described by Snell's Law. Examples of 
some simple ray paths are shown in Figure 2.1. A gene- 
ralyvsed raysinciudes the head waves and interface waves 
associated with such a geometric body ray. The response 
from each ray is found by solving a wave equation of the 
EOmine 246) 9200) OF (2.0) for cach ray esegqmene. on ule 
generalised ray, with the appropriate potential and 
velocity for that layer and segment wave type. The dis- 
placement and traction boundary conditions are then a 
applied at each interface in the order encountered by the 
bay patn, to obtain the particular solution for tie Speci— 
fied ray. 

When a P or SV Segment of a generalised ray is inci- 
dent on an interface in a solid, its energy is divided 


among the four possible reflected and transmitted P and 
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oVi ray segments shown in Figure 2.2. In a fluia, the Ss 
waves are absent. An incident SH ray segment will gen- 
erate only the reflected and transmitted SH waves. The 
energy going into rays of type different from the next 
segment of the given generalised ray is removed from 
the ray at each interface along the ray path by the 
boundary condition solution. This energy is of course 
included in other rays of the ray expansion. The gen- 
eralised ray expansion in a model of K plane homogeneous 
elastic layers consists of the infinite number of ray 
paths consistent with Snell's Law. The impulse response 
SOLUtELON tO the equation of motdonm v2. 1)%. found from “a 
generalised ray expansion, is a series of infinite spikes. 
When convolved with a pulse, obtained from the source and 
receiver parameters, the.solution is equivalent to the 
complete solution obtained by any other method (e.g. 
Haskell matrices). In practice only a limited number of 
rays can be calculated. The generalised ray sum must be 
ordered so that the impulse response converges when con- 
sidered as a generalised function or distribution 
(Pepoutic.) 1902, 90.270), Only. therminimum of rays of 
sufficient amplitude to construct the final seismogram to 
within an acceptable amplitude error are then selected. 

When the model parameters vary slowly with depth, 
rays could be selected in order of increasing number of 
reflections, since the small reflection coefficients are 


the dominant factor in determining ray amplitude. Since 
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the amplitude of a ray decreases at every interaction 
with an interface, another procedure is to select rays 
in order of increasing number of ray segments. This 
latter ray selection process can be done systematically 
using kinematic and dynamic group classification (Hron ,; 
1971). Any two generalised rays having the same number 
of path segments 2ny in each layer k of a model of K 
layers have identical temporal properties. For rays 
with no conversion of energy between P and S propagation 
modes, the number of such kinematic analogues is (Hron, 
ae Se del) 


K 
N i oe : C27 lela) 


Sabena gion oer 
GN ai ee eet Nya] 
ae is the binomial coefficient (Abramowitz and Stegun, 


ee Oe oe 2.) 


o = £!/m!(L-m) ! . (254-2) 


These kinematic groups can be subdivided into 
dynamic groups, all rays of which have identical ampli- 
tude properties. All these dynamic analogue rays have 
the same number My of reflections at each interface from 
incident ray segments in layer k. The number of such 
dynamic analogues in each dynamic group is (Hron, 1971) 
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The summation of the ray expansion can then be 


assumed to converge as a distribution when the rays are 
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arranged in dynamic and kinematic groups, and the 
kinematic groups are summed in order of increasing 
number of segments per ray. 

There are two advantages to the Hron ray classi- 
fication scheme. First, the ray expansion has been 
shown to be complete when no P-S wave energy conversion 
occurs (Hron, 1971). All important rays can be easily 
included automatically in any ray-method solution. 

The second advantage is the saving in computing time 
due to grouping of properties. Lea sOlweLOneOL= Lic 
boundary conditions is done only once for each dynamic 
group. The temporal properties are determined only 


Once Lor=all N rays (2.11) in each kinematic group. 


KIN 


a3¢°5 


re eine 


7“ 
peahpirnnVotiCa yy oe 

‘qhivseriee Aeneke whew 30-00 i ae Ja ae 
¢itese sd’ nss ayes shedzoqm? ITA «(eva amis exuwed 
inoisufoe boilsam-ysx ye ni yilsoktemasin babu font 
emis paitugmos mi privsa siz Bi- - apsdiisvbe bnovse ost 
BY Go noiswloe SH? .zatsaSqoxc 35 eitapoze oF sub 
oimenys ddae to Soro ylne snob 2i anoidibros esnturbeddl 
yino beaimetsh sxe asidasqgoxq fsiogme+ eal quote 
amoxze oigsmetix dosevai (11.8) eysa »..¥ Lis set sage 


25 


CHAPTER 3 


THE CAGNIARD-DE HOOP METHOD 


3.1 The Transformed Wave Equations 


bet. xX and fF be @ pair of potentials in (2.5) 
satisfying a wave equation of the form (3.1) ina 
homogeneous medium 
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The problem is to find the response ,(t,r,6,z). Follow- 


ing. Cagnvard’ symethod: (Cagni ard, i119 39) 1,0 thersolution as 
found by transforming the differential equation varia- 
blest fir sland: 60‘solvitig Hori the transform y(S;p,2%,2), 


then taking inverse transforms to obtain y(t,4r,0, Z). 


S.teL The Source Function 


A realistic seismic source term is a point dis- 
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pace eee ae 2 a (352) 


The factor 1/r is necessary to make the total source 
strength finite when integrated over volume. The form 
of the force f is found by substituting (302) into the 
force equation (2.5). The differential operators V 


and yx in (2.5) operate on the field point and, source 
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separation (r-r_), where the variable is the source 


is 
=O 
location xo: (OA prime represents a derivative with 
respect to a source co-ordinate e.g. 6'(r-r.) = 
dé(r-r.) /dr.. mnie vector Feis) (0,0,\%):. 
Pie .LOLcesdistribution. from, lis chen 


6(r-r) 5 (6-6) 6 (2z-h.) 


£,= (At2u) VP = -2 (A+2u) Fo (t) | ———+__"___-_ 
6 (r-r,) 6" (8-8.) 6 (z-h.) 6 (r-r.) 6 (0-6) 6'(z-h) 
CS ae) ! xr - 


The delta derivatives then give force couples directed 
radially outward as shown in Figure 3.la. This type of 
explosive source generates only P waves. The force dis- 


tribution associated with M,, (2 Sass 


iu _ rrr) d ie 
fe) 
Os Zit F(t) ‘ 
In this case the delta derivatives give a horizontal 
double rotational force couple about the z axis (Figure 
3.1b) i? / This ‘can-cause.only SH iwavesi: 
The force distribution arising from M, (253 2ouS 
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b) HORIZONTAL SHEAR (SH) SQURCE 


Figure 3.1 


The Three Source Force Distributions 
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h) TOTAL SV FORCE 


VERTICAL SHEAR (SV) SOURCE 


Figure: ssl) (continued) 
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The force couples due to each component of (3.1.4) are 
Shown an higures.%.iL¢ to 3.1lg, and the total torce 
distribution is represented in 3.lh. The force is 
always in a vertical plane, and generates SV waves. 
Other types of sources having azimuthal varia- 
tions, or having finite dimensions were not considered 


in this study, although the extension is not difficult. 


3.1.2 The Transformations 


The Onagin Ls located Lor Convenience directly 
below the source on the bottom interface of the layer 


Ko containing, Che Source, 


First apply the Laplace transform of time t to 


both sides of the wave equation (3.1) 
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This is followed by the finite Fourier transform of angle 
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and the Laplace-Bessel transform of horizontal distance 
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The transformed wave equation is then given by 
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The velocity v may be either q (2.9) or g (2.10). The 
vertical wavenumber of Ny 4k in layer k of the model and 


mal 
for=ene= 4 : ray segment of a generalised ray is 


2 HE Z 
vik ne 2 Pral* (ane 
k 


3.2 Transformed Response for a Generalised Ray 


For each ray in the expansion, the transformed 
wave equation (3.6) must be solved in each layer of the 
model. In each layer k, the transformed solution of the 


homogeneous form of (3.6) is 
X(S,;Prh,Z) = AL exp (“Sn,5,2) An exp (Sn.,32) . (3.8) 


The constants Ay andpA, contain Lhe, amplitudes: of the 
waves travelling in the positive and negative z direc- 
tion in layer k. In the layer Ko containing the source 
fe, at) depth Ao radiating energy of wave type eee the 
solution must include the additional particular solution 


found by the Wronskian method (Morse and Feschbach,1953,p.530) 
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This solution represents energy radiating from the source 
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The amplitudes A, and A_ are found for each ray 
segment by satisfying the boundary conditions at each 
interface in turn, starting from the known source radia- 
C1onicerm (S.0). 

For the general P-SV reflection-transmission 
problem at interface k, there are four equations of 
the form (3.8) describing P and SV waves on both sides 
of the interface. The 5H problem is much simpler so 
will not be done here. For simplicity the layers are 


labelleds lWeands2.. in layer 1: 
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Similar equations apply in layer 2. The eight factors 
Bee and 8) sag are the wave dmplitudes’ to be deter— 
fs) (2) i) (4) 
mined. The other constants are introduced so that a 
unit wave amplitude represents unit energy flux in the 
ZeGLreCciuon. 

At the interface, the displacement and traction 


are continuous. These quantities are conveniently 
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oF ses 
aV. 
Toe (3.11) 
ie rats. 
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ox Pore ‘ 
u; and Oa uere displacement and stress. This vector V 


is) found by applying the transforms (322) (3.4) and 
(3.5) to the displacement equation (2.5) and to the 
constitutive equation (2.2), and substituting the trans- 
formed potentials (3.10). 

In layer 1, the transformed stress displacement 
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In layer 2 V, is given by the analogous equation. The 


Z 
origin has been shifted vertically in (3.12) so that 


the new z coordinate z' is zero at the interface. 
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are displacement-stress eigenvectors for a wave of type 
P (S) of unit vertical energy flux travelling in the + 


(=) 2" cirection in layer 1 (2). 
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The stress-traction boundary conditions imply 
V, = V, at the interface. The amplitude P, | or Ss i 
of the incident ray segment is known Pee ike source 
term at the first interface, or iteratively from the 
previous interface boundary conditions. The amplitudes 
of the other three possible incident waves are taken as 
zero. The four remaining amplitudes of P, j] and Say 

(a2) (=) (2) 

are then the amplitudes of the four waves generated 
at the interface, as shown in Figure 2.2) 

The) ratio of the amplitude of thej reflected or 
transmitted wave % to the incident wave at interface k 


at the end of ray segment j is called the reflection 


transmission coefficient KS om (I), Expressions for these 
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coeftiicvents areygiven in Chapter 5. 
Phe derinition of possible values of § and! m 
below is similar to notation used by Cerveny and 


Ravindra (1971) 


A P wave above interface 
S wave above interface 
3 P wave below interface 
4 S wave below interface 


The one wave amplitude corresponding to the next 
segment j+l of: the generalised ray is then taken as 
the incident wave amplitude in (3.12) for the next 
interiace). ~The origin isvagain shifted so, that 2. is 
zero at the next interface, and the boundary conditions 
again applied. 

At the receiver in layer k, at depth h, the last 
segment J has wave type Vie The solution £or the géen— 
eralised ray is 
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C(p) is the amplitude change from source to receiver to 
satisfy the boundary conditions i.e. due to reflections 


and transmission at interfaces, divided by Ny lk ° 
el ke) 
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Oyis the phasesftactor introduced by origin shifts ito 

Make 2 Zero au cach anterface in. turn.) This represents 

the total phase change of the wave along the propaga- 

tion path of J segments 

J 

2S A Nik 4k ; (2216) 
j=l 

The factor Osx is the vertical distance travelled in 

layer k by the ae ray segment. At the source, diy is 


Ene source. cepca Ay if the first segment goes downward, 


and ike = ho) if the segment goes upward. Here se 


is the thickness of the layer containing the source. From 


‘ The other 


the final interface, at depth ny dees een 
Boy, are the thickness of layer k. 


oo The Complete i rans tormea Ray So ittrion 


The equation (3.14) 2s the transtormmed response 
for a single generalised ray. The complete ay expan- 
sion for P waves in a modél of K layers is given: by 


(Chapman, personal communication) 


far .&) 


a2 ime Tawi ‘ts iii el 4 


ve a” (Figgne se I 


-buarevct, w@e Wetese 2)" wos °f a) Pring 
ar 4t0% ,.cLsaq a - i ea » eon). =_/! = 


-SOFa to) wis er ee 


wy =s 


Fos" 
saGsn wi? 


s 
. ne 


verse mad) f 


wera Amd yt ere= S 
nbere Gee Oooh Gese per 


aoe 


Cass ¢! somiretn 
negagnn G4) winks view ws 1 wen =ae Tuts 


- Sreeerr fe: 


a Lo 
= Sept ay ~-s @ 
mele a wee <a (Dw) 
Cl ald Vet oem |) 


bon am, inn? 


2 


eihesje= < Si dean! ant? 
i jeu" 


murat: aff 


ri AY Gis? . ha 


g 


\e 


7 
es qt é¢ weyet 


= 


—pese old 


: + 
ut Sih 
& = 
ogup Vi's eri uc aly eno ae: 
=> 
2 8) Shae Fa ren a2 i 
ed 
4 : & SH) 7 ab wie PT 


5 
ee ee a: 
age wil? : 


fi €HGA_8 ae 


38 


y ae 
e Bale) es K 
s ee pa 
MUSZOr 2) 2 S| L exp as are ey 
orto 
il. ==0) 
n,.=0 
n,-1 
Kal nea 1 n 
‘i (,833) ; ce la 
m,=max (0,n,-n.) 
m,=max (0 ,n,_4-Ny) 
K-13 ns rip m. n.-m. me -n.t+m. 
Ek i+l Be i bg te ye tl alia ia 
ee ene ae ea ay (533? 
11 Bf oh) ee 
Caer) 


The numbers ne, mM; and Coa are defined in Chapter 2, and 


the are reflection-transmission coefficients. The 


k> om 
outer summation ranges over kinematic groups whose 
temporal properties are determined by the exponential 
factor. The inner summation ranges over the dynamic 
groups within each kinematic group. The amplitude pro- 


perties of each dynamic group are determined by the 


reflection and transmission coefficients ke om* 
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The ray expansion is more difficult when energy 
is converted between P and S modes, but rays with one 
or more converted ray segment can be included in the 
automatic computation of the ray expansion (Hron, 
personal communication). In many studies, rays with 
more than one converted segment have arrival times 


outside the time window of interest. 


3.4 de Hoop Inverse Transforms 


The inversion of the transformed solution must 
bemdoneeforsone raycofeeach dynamicigroup. in? thetray 
expansion A6t4a7) 

AfEterstakingrthevinverser transforms (8.5): Of «p 


and (3.4) of j}2, the solution (3.14) becomes 
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The time domain solution could now be found by taking 
tne inverse transtorm (3.4) Of Ss, and evaluating the 
double integral to obtain the response at each time 
and position required. However, it is’more efficient 
to follow the method introduced by de Hoop (1960) and 
used by Helmberger (1968) to avoid both integrals. A 
Change of variable from p to £ is used to put the right 


Nandesice of (3.18) into thertorm (3.3). 0f the forward 
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Laplace transform of time. The response y(t,r,0,z) is 
then extracted from inside the integral by inspection. 
First the asymptotic form of the modified Bessel 


funectvon 1s substituted into (3.18) 


| TT 
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This approximation is valid for early motion teleseismic 


studies because frequency s/i and range r are large. 
The approximation is not strictly necessary to obtain 
the solution, but greatly simplifies the procedure. The 
complete result is obtained by the more difficult method 
of double integration mentioned above. 

In order to remove all dependence on s in (3.18) 
other than in the exponent, assume a source function of 
the form 


F_(s) = F_/vs 
2 —) (3.20) 
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Migsechoice of source Eunction 1s required by the nature 
of the po@nttsource problempebut’ caneeasily be removed 
and replaced by any arbitrary source by a convolution 
after the total response x(t,r,98,z2) has been found. 
Thencsubstitubaonse(sttojeand (3.20) and the change of 
variable 
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put (3.18) into a form superficially resembling a 


forward Laplace transform of t 


F p(t)=i° 
ee roy cep mie io dp 
4 r Zz ES (C(p) vp Fn) EXP ( st)dt 


1 p(t)=-ie 
(S22) 


Since the Laplace transform of a: physical quantity 
Vie, G,0,2Z) must be real, 


(C(p) Vp SP) exp (-st) dt 


Sahat ok (3.23) 


ele Oak = 
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using the principle of reflection (Titchmarsh, 1939). 
If (3.23) is to be treated as a Laplace transform, the 


variable t must be real. The contour defined by 
p = p(t) Im(e)e =20 (3.24) 


found by solving (3.21) for p, is called the de Hoop 
contour. The de Hoop contour cannot be written expli- 
Citly in a multilayered model, but usveasily located 
numerically z 

When the p plane integration contour of (3.23) 
from 0 to tio can be deformed into the de Hoop contour, 
@3i2 3)? ciisHia'tva ltd ‘Laplacelttransfiorm,; and? y(t pr (o-p2 ees 
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The main features of the de Hoop contour are shown in 
Figure (3.2). he contour can be started at t = 0 on 
the negative real p “axis, ‘and ol lows ‘the realetmaxis. 


The response is zero here because the function 
a 
(C(p) ¥p ge) 
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At p = 0, (3.26) is positive. However, at p= Ly eae 
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dt/dp goes to -~, since all nik are real and positive. 
Necessarily at some point in between, called Pray! 
dt/dp =0. This means that the complex function t(p) 
always has a saddle point on the positive real p axis. 
The orientation of this saddle point is given by the 


second derivative 


2 J 


Gat Pha ies: 

—z~=- } d.,/ve no. ‘ (3.28) 
dp? vey) af eee ae) ame Viayite 

At Pray , the second derivative (3.28) is real 


and negative, so the de Hoop contour, being a steepest 
descent path, enters the first quadrant at 71/2. At large 
values of p the de Hoop contour becomes asymptotic to 


the line 
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Arg(p) = 1/2 ve 
rg(p) = qT = arctan |. rele : (eo) 
spa 


The layer index is k and the ray segment number is j. 

the gen tegration: path in p(3.23) can be distorted 
to, the del hoop contour! The antegwal is Zero, onthe 
arc at infinity joining the imaginary axis to the de 
Hoop contour, and no singularities of the integrand are 
crossed in the contour distortion process. 

The seismic arrivals are correlated with singu- 
larities or the solution (3.25). The infinite discon 
EinuLty Crom dp/dt tat the saddle point is tne rerlected 
geometric arrival. In the vicinity of the saddle point 


Pray’ the response (3.25) is approximately 
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H(t) is the Heaviside step function: H(t)=0 
H(t)=1 


Head waves arise along branch cuts of "vik from 
tomcOuDranch posits p=1/v, on the real p axis to the 
Vett of the saddle point, for any layer k touched by the 
ray. Along the branch lines, C(p) has an imaginary part, 
So the sesponce (3.25) Us, monzero. » The first motion from 


(3.25) at a head wave is approximately 
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X(t | + 6t) - aay = Fy (2tr) K, (dt) * 
(33315) 
y 3 - 
K, = lp Boy See sees) 
(Sp) P=Pp 
1 


Interface waves are associated with poles of 
C(p) near the de Hoop contour. This de Hoop contour 
is found once for each kinematic group chosen from the 
Lay expansion’ (3717) The *solutron’(3e25) is then 
evaluated along this contour once for each dynamic 


group in the kinematic group. 
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CHAPTER 4 


THE EARTH FLATTENING TRANSFORMATION 


One purpose of this study was to examine seismic 
Signals near the shadow edge 8 of the earth's eae 
core. The three imayor types oF* signal of iampomtance 
at tiese Yanges are 11lustratedeby rays in Figure: 4.1. 
In the illuminated zone, 6< oF the direct wave 1 and 
the core reflection 2 exist. Beyond the shadow edge 
GS : there is just one diffracted wave 3. 

The Core shadow edge™is* observed at ee 265 
(Sacks, LOOG rs At this large range, the curvature 
of spherical interfaces cannot be neglected. The 
Cagniard-de Hoop solution cannot bevapplied directly. 

Gilbert and Helmberger (1972) have shown that 
the Cagniard-de Hoop method can be applied approximately 
to a spherical model by using spherical coordinates 
centred at the sphere's centre. After approximating 


ModLiaed Bessel functrens by (HErdelyi set al, 1953), 
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they obtain a solution of the same form as the plane 
Payer solution (3.25) , multiplied by a geometrical 
Spreading factor (0/sin 8), Since p is the transform 
of an angle 6 instead of a distance, the plane wave- 
numbers (3.7) in C(p) are replaced by spherical wave- 


numbers 


f 4 She op Cor) 
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k 

This method of Gilbert and Helmberger can only 
be’ used for reflected “rays (2°in Figure 4.1). For ray 
segments of the direct wave 1 travelling nearly horizon- 
tally, the Bessel function order number } is small, and 
the approximation (4.1) is invalid. Although a solution 
of the equation of motion exists it cannot be found by 
the Cagniard-de Hoop formulation. 

The diffracted wave 3 also cannot be obtained by 
the Cagniard-de Hoop method. Beyond the shadow angle 
Our the saddle point equation dt/dp = 0 has no solution. 
The saddle point and a de Hoop contour do not exist. 
Knopoff and Gilbert (1961) have found the asymptotic 
solution deep into the core shadow by directly invert- 
ing the transformed solution analogous to (3.18), but 
wave methods must be used to investigate diffracted 
Signals near the shadow boundary. 

It is known that direct and diffracted waves 


exist. The Cagniard-de Hoop method is just unsuitable 
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for the problem. The approach used in this study was 
to apply on earth-flattening transformation to the 
model parameters, changing the spherical surfaces to 
plane surfaces, in such a way that the kinematic response 
(arrival times) due to an impulse disturbance was 
unaltered, and the dynamic response (e.g. displacement 
amplitudes) was minimally affected. It then followed 
that the response of this flat model was identical to 
the response of the original spherical model, to within 
the amplitude error associated with the transformation. 
The flattened model was then approximated by a stack of 
plane homogeneous layers. The complete response of the 
flat model at all ranges less than TX, could then be 
found by the Cagniard-de Hoop method of Chapter 3, merely 
by including the spreading factor (@/sin 6), Chts 
response of the flat model was also the response of the 
original spherical model of radius ro: 

This method has been used by Muller (1970) and 
by Helmberger and Wiggins (1971). The optimum form of 
power law earth-flattening transformation has been 
determined by Chapman (1973, in press). He found that 
the transformation used by Helmberger (1972) is optimum 
in a fluid, wnile the transformation of Biswas and Knopoff 
(1970), shown in Figure 4.2, is optimum for SH, SV, and 
P body waves in a solid. This transformation was used 


in this study. Quantities on the left are flat model 
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parameters 5* *Those’on ‘the *right Vare “from "the @spherical 
model. x and Z are horizontal and vertical distances, 
and V and RHO are velocity and density respectively. 

The geometric effect *of *the veloerty “transiorma- 
tion is to increase the ray curvature, compensating for 
the decreased curvature of the model surfaces, so that 
the arrival times are not altered. This is illustrated 
ny Figure 4.3. 

Chapman (1973) has written the transformed wave 
equations and continuity conditions for the spherical 


model in the form 


Ve is the transformed spherical displacement-traction 
vWeCLOl,adnaLogous €On(3.L1):, W is the transform of the 
source terms, and M contains model parameters and boun- 
dary conditions. The earth-flattening transformation T 
is applied to M and written as an expansion in (ashes 


Si is the time transtorm of (3.3), and a is the core 


radius 


(0) 


M is the matrix used in the flat model formu- 
lation. The approximation in the transformation T is 


measured by the higher terms. 
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the tCransrormation in BPigqure 4.2 1s exact, for 
HOvVe waves. CLLOrs OF 0((xps)~7) ake Present ano 
and SV amplitudes, and of 0((uPps)~*) in P amplitudes. 
First motions are unaltered by the transformation. 
Some pulse distortion may be expected due to the fre- 
quency dependence of the error terms, but the effect 
is negligible for body wave studies where frequencies 


are high. 
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CHAPTER 5 


NUMERICAL CALCULATION OF SEISMOGRAMS 


Side Outline 


The calculation of a numerical seismogram is 
started by generating an earth model of plane homo- 
geneous layers. This may be done by defining the 
flat layers directly, or by applying the earth- 
EVaLtcening cranstormation (Chapter 4) toa» spherical 
model, then approximating the resulting flat model 
by plane homogeneous layers. It is also necessary 
to specify the source and receiver characteristics, 
the time window ae to oe ey, the digitizing interval 
DT, the number of rays in the ray series, and the 
accuracy required in calculating the rays. 

Kinematic groups, defined by sequences of inte- 
gers n, (Chapter 2) and dynamic groups, defined by 
sequences of integers m, are generated automatically, 
in order of increasing number of ray segments, up to 
the preset limit. 

The kinematic properties of each kinematic group 
are found by locating the saddle point and branch points 
of the response (3.25). If the response at these impor- 
tant times is greater than a preset value, the response 


is found at all times within the specified window, by a 
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Mimi tegsanpeergemethnod. Starting: at the first arrival, 
steps are taken along the p-plane de Hoop contour, and 
the total group response is found at each point. The 
stepisizesis controlled by .dy/dtand d*y/at* so that 
small steps are taken near large variations in response, 
and larger steps elsewhere. .The reflection coefficients 
kom are obtained at each interface for each de Hoop 
GContourm point in. onder sto, .evaluate .the iC(p).in (3.25). 

The procedure is stopped when t(p) is beyond the 
end tase of the preselected time window, or when an 
estimate of spr indicates that the Bessel function 
approximation (3.19) is no longer valid. 

The response of the kinematic group is then 
smoothed and digitized at the desired time interval 
by ajconvolution with a triangular window. .This result 
is added to the total response from all previous groups 
in the ray expansion. 

The spectrum of the final ray sum is obtained by 
aulact sROURLOE,transtorme (hoe .T.) (Cooley and Tukey, 
L965 ).,enultiplied by (s)? =- LW) <.to,deconvolve the 
SOULCeeSUNCLION s3.20)i,mandumultiplied by the bbl. 
of the source and receiver. The inverse F.F.T. then 
produces the desired numerical seismogram. 


This seismogram calculation procedure has been 


arranged by C.H. Chapman. 
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5.2 Barth Models 


The models used consist of a stack of up to one 
hundred plane homogeneous layers on a half space. 
Longitudinal and shear velocities, density and thick- 
ness define each layer. Since energy restricted to 
surface regions often arrives later than the time 
window of interest, up to fifty layers at the top of 
the stack can be specified as a transmission - only 
section, represented by DMODEL and UMODEL in Figure 
5.1, i.e. the ray summation does not include 
rays with reflections in this section. The section 
merely creates the correct geometry and amplitude for 
deeper-penetrating rays. 

The structure in this transmission section may 
be different under the source and the receiver. This 
is useful for paths from oceanic earthquakes to con- 
tinental seismographs. The solution for such a non- 
symmetric model cannot be found eeacel yi tnecretsealiys 
but when no head waves or near-horizontal rays are 
generated in the transmission section, the response 
obtained by the method of Chapter 3 is adequate for 
Practical purposes. 

A ray expansion is then generated automatically 


in the continuous layers labelled MODEL in Figure 5.1. 
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5s SueAUtomacic Ray Generation 


Before a ray expansion can be generated, the 
source location and wave type P, SV or SH, the maximum 
number eed of segments in a ray, and the maximum pie 
of ray segments of converted wave type (P-SV) must be 
specified. 

A ray with no conversions in a model of K layers 


is defined by the series of integers n, giving the 


k. 
half number of segments in each layer up to K. The 
kinematic groups are generated by incrementing ny in 


the shallowest possible layer k to obtain a new sequence 


such that 


The dynamic groups within each kinematic group 
can be found by incrementing the numbers my of reflec- 
tions at interface k from incident segments in layer k 
to obtain new ray paths consistent with the set of Nys 
The number of analogues in each dynamic group is then 
FOUNGTELoms (2.13) < 

If rays with one converted segment are also to 
be included in the ray expansion, then each segment in 
turn of each dynamic group is replaced by the segment 
of different type, to generate a new dynamic group with 


one converted segment. More complicated algorithms can 
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be used to generate groups with more than one converted 
segment (Hron, personal communication). 

tierray rerlected trom the top. of the rerlecting 
zone (MODEL of figure 5.1) is also included when the 
transmission layers exist. 

Rays going upward from the source can be inclu- 
ded when the source is in the reflecting zone MODEL. 

When the response for refracted waves in a nearly 
homogeneous section of the model is required, (e.g. P 
phase at ranges less than 100°) a practical ray expan- 
Sion can be generated to include only once-reflected or 
primary rays instead of the dynamic groups described 


above. 


2.4, ce de, Hoop, Contour 


Andel Hoophcontoursinvthe p plane (figure s3.2) 
must be located once for each kinematic group of rays. 
The saddle point p and branch points are located 

ray Ph 


tLset.e AnVupperaiimi toP2eon Pray is given” by Sh a 


t/a tS mne(em27).athestvetocities Vv, in all layers 


touched by the ray are then tested for head wave branch 
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tnen sche branch point 1s tothe left of the saddle 
and a head wave exists. A lower limit Pl for the 
saddle point is given by the maximum of all P, = 1/v, 


Saeisty ing 205 i jee o ke 
fdt/7a <0 
7. Plo=p, (5.2) 


then this branch point is to the right of the saddle 
and no head wave exists. An upper limit on the saddls 


point is given by the minimum of p and all 1/v, 


max 
Saeustyving (5.2). 

The saddle point Pray is then found by Newton- 
Raphson iteration (Abramowitz and Stegun, 1965, 25,2.28) 


Frome ue lore bP 2) 72. COntour Steps oF 


Sp =-at/ap /a*t/ap? (5.3) 


are taken until 6p is less than a specified error limit. 
Expressions for the derivatives in (5.3) are given by 
(Sie26 anc. 2c). « 

The saddle point value Pray is, 1dentueal to the 
standard ray parameter p as used by Bullen (Bullen, 1963, 
p.109) in geometrical ray theory. 

The importance of this kinematic group is then 
estimated by evaluating numerically the amplitude change 
(3.31) at each head wave Phy with a small step of ép 


along the de Hoop contour, and evaluating the amplitude 


near the reflection using the approximate response (3.30). 


(8.0) 


elbise si3 30 sipin saz ‘cot ak vaiog donaxd eindt ats 
efbbss oft no simti rage cA .ateaxe =vsw basa on baa 
gut ifs one , 2g 6 sacl Sit ye mevip at saiog 
(8.2) oaiytedtee | 
“foiwek ya Basel ‘nenz et = g tntoq slbbse ont 
(88.8.8 Jeaed ,mupess bas yell molsjstagt noadget 
26 eqatze suotnoD §=.S\ (89 + ES) moat, 


te} Sap\ 25 \gp\sp+= qa . 


-timil tozuxe Ssittooye s nsdg eesk eb ge lishe wedtes e158 
yd davyip sas (£.%) ab asvisevixeb sas 102 noo tees3qxa 
(BS 48) Bas (884) 
yas? Ssulsy gmiog sibbee edt | li 
E3EL \MeliUs) Asifed yd peau esq wesomexsg Yor Exebanse 
syaasit yea Lankenres mi (ROL. “a 

mois <i qlorp SitementA eine to onic 7 oT, 
epaside eee Lae oat yiisodramsa | ist 
gi..26 gese Ifsne 6 cat fc pe os nial ie. E) 
Sbugitqna: ods prbtasiate b tre 1088 sab mae eds patie. 


ioe tesa stan ates 


6H? cd Lsoitnebi al 


i. > 


— 


60 


Pulses with dominant frequencies near the peak 
frequency of a combined source-receiver spectrum are 
expected to be important in the final seismograms. 

The early motion amplitudes above can be compared in 
an approximate way in the frequency domain at the peak 
source-receiver frequency Te ee Since the head waves 


1, 
and the reflection are approximately K,6t* (Bias) vane 


-t1 
K,6t ° (3.30), their Fourier transforms are approximately 


Lit =, 
proportional to 2K 5 evan Kyu phe Fie K, or any of 


K5/2W ax is greater than a specified threshold amplitude, 

evaluation of this kinematic group response is continued. 
The response calculations (3.25) are started at 

the first head wave. At each head wave, the response 

at Ph. and Ph, + oP is already known. The subsequent 

points p(t) used on the contour can be chosen by using 


a quadratic representation of X(t,r,6,2) 


a ax a i ;ax te 
5 ee 
J 
(5.4) 

The derivatives dxy/dt and ay /at? are approximated 
by finite difference expressions from the previous points. 
Since three points are necessary for a finite difference 
estimate of d*y/dt’, the second point after a saddle or 

2 


branch point is estimated assuming d*y/dt = 0. A maximum 


acceptable fractional amplitude change C between adjacent 
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points then determines the largest allowable time step 


Gt = 3) in (5.4). The p-step on the de Hoop con- 


see 
tour is then 

(pe =n) oe ee Coe = t2 - 

Deemer) te bv Se cece (oy) 

Beyond the saddle point Pray! the de Hoop contour 
is complex. The time Sirol for the next point is estimated 
as before, but a Newton-Raphson iteration process must be 


= 0.) By idetingtion p. 


applied. tO. pp. to make Im(t. j+1 


+1 4+1) 
then lies on the de Hoop contour. 
When each point Paad has been located on the con- 
tour, the response (3.25) is evaluated at Pay? multi- 
Dited by the dynamic group multiplicity (2.1393. and 
summed over dynamic groups to obtain the exact total 
kinematic group response at t = je ail Time ot is 
determined exactly from (3.21) after Sn has been found. 
Point selection is terminated when SiPay a! is later 
than the end Sa of the specified time window. Since 
the Bessel function approximation (3.19) is invalid 


for small arguments, an estimation of spr should be 


used to determine the validity range of (3.25).° If 


Poa EONS S| - tray? =sspre~ L0ass, (5.6) 


then point selection is terminated, and the amplitude of 


the ray response is tapered smoothly to zero at the edge 


| = a4 eee 
(2.2) : Pd aa 7 Ine "US - re? 
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Ofnche gregionsoravalidi ty: 

This method of point selection results ina 
high point density near major arrivals when the 
response varies rapidly, but a lower point density 
elsewhere when the response varies slowly. This 
enables accurate seismograms to be obtained with a 


Minimum of calculation. 


5.5 Reflection Coefficients 


At each point Py on the de” Hoop contour, all 


the reflection-transmission coefficients kom (Chapter 


3) needed for the kinematic group of rays, must be 


calculated. The boundary condition equation Vy = V5 


FLOM Sel 2 tp ls Cs 


+ P_lV Pas acme WY a 


Fea EDT eer 1 461 ee ee ae 


pay as fs (5.7) 


=o pp Oconee ea 


can be solved analytically to obtain general expressions 


Por Sat an arbitrary iImterface k. These expressions 


k> om 
are shown in Table 5.2! for P and SV Waves at a general 
solid-solid interface. The interface index k has been 
Omlecea for Simplicity. 

The terms in P_j, S_47 Pio and S,., si ed ee 


represent waves incident on the interface. When one Ot 


these four amplitudes is taken to be unity, and the 
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other three to be zero, the four remaining amplitudes 
Ploy Siar Plo: and S_> are the corresponding reflection 


and transmission coefficients S The set of sixteen 


&m* 


coefficients can be found by a matrix method (Dampney, 


LOA. 

2 Sia Se ee) ee eee (5.8) 

= —+Pl—+Sl1—-P2—-S2 —-Pl—-S1—+P2—+S2° ° : 

The elements of the four by four matrix Sere Som 

and each V an ee iS a four by One columnivoL 4a Lour 
C=) (Sy 20) 


by four matrix, represented by brackets. The expressions 


inelable<s. fLores m are obtained by applying Cramer's Rule 


x 
(Nomzzit, £960, 0p. 179) to (528)..—" Because of the normalzza— 
tion Usecde nels. 10); thesmatnix S is isymmetric so that 


S.. = S.. and only ten independent coefficients need be 


The coefficients for interfaces involving fluids 
can be obtained by letting the appropriate u and 8 in 
Table 5.1 approach zero, provided denominator factors 
of Ne have been cancelled by rearrangement of the ex- 
pressions, or by using fluid boundary conditions in 
Chapter 3, i.e. there is no shear stress in the solid, 
and tangential slipping at the interface is permissible. 


Normal displacement and traction are continuous. 
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the’ coerficients for the simpler SH reflection 
problem are™shown in Tabre* 572. 

When the model contains many layers, the trans- 
mission coefficients at near-surface interfaces vary 
slowly with p for deeply-turning rays. Great computa- 
tional savings may be obtained with an acceptable small 
amplitude error by calculating the coefficients for the 
top layers once e.g. at the saddle point, and using 
these, Cocttrcients at all” points on the de Hoop contour. 
the coefficients at all interfaces within a specified 
number of layers of the rays turning point should be 


recalcudated at each point on the contour, however. 


5.6 The Smoothed Time Series 


The response as calculated above at variable time 
intervals\.is difficult) to manipulate for the addition 
Of Utne ray Expansion, sor for convolutions or fast Fourier 
transforms. It is more convenient to have a time series 
sampled-at-unitorm=time-intervalseDT.--in addition, the 
response found in Section 5.5, although exact, contains 
high frequency components that would be aliased by a 
fast Fourier transtomm (Nyquist, 1928). “A triangular 
window w(t) (Bartlett, 1950) (figure 5-2) 1s convolved 


with the response from Section 5.5, and the result 


digitized at intervals of DT. When the window width 
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WNDW = 2DT, the first zero of the window spectrum 
(figure 5.2) coincides with the Nyquist frequency 
Tt/DT of the digitized seismogram. This means that 
frequencies greater than 1/DT have been effectively 
removed from the seimogram and serious aliasing is 
avoided. 

The smoothed, digitized response for the 
kinematic group is then easily added to the total 
time series sum from previous groups in the ray 
expansion. This total sum is easily manipulated 


WEtEhnstchne Last Fourier transform, (F.F ol.) 


Pinie walt Arbitrary Source and Receiver 


The Cagniard-de Hoop response (3.25) has assumed 
aysournce function F(t)” = Be (3.20). Glt tsanecessary: 
to deconvolve this source and convolve the result with 
a realistic source mechanism, and a receiver response 
function, to obtain a realistic seismogram. These con- 
volutions are done by means of fast Fourier transforms 


(PoP t.)s and multiplication an the Grequency domain: 


Multiplying the F.F.T. of the seismogram by 


1 ail 
=a 7% 
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(iw) * removes the §S = (iw) spectrum of the source 
(3.20). The source spectrum of an earthquake or explo- 


sion can be obtained by a deconvolution using the 
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response at some seismic station where the seismogram 
is simple. Multiplication by this source spectrum and 
by the receiver transfer function then produces the 
final %seismogram spectrums, An inverse F.F.T.acthen pro- 


duces the final time domain numerical seismogram. 
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CHAPTER 6 


NUMERICAL SEISMOGRAMS IN A HOMOGENEOUS EARTH MODEL 


FROM THE ILLUMINATED REGION TO THE DEEP CORE SHADOW 


G.2 l£ntroctetion 


Numerical seismograms at all regions of a 
homogeneous earth from the illuminated zone across 
the core shadow edge and deep into the shadow have 
been obtained by a wave solution method (Chapman and 
Phinney, 1972). Gilbert and Helmberger (1972) dis- 
cussed an approximate form of the Cagniard-de Hoop ray 
method that is valid for computing the wave reflected 
from the core in the illuminated region. Knopoff and 
Gilbert (1961) obtained the asymptotic form of the 
Laplace-transformed diffracted pulse deep in the shadow, 
in terms of Airy functions. By then performing inverse 
transtormations of the form (3.3) .and (3.4) ,. they found 
the time domain pulse deep in the core shadow. 

The signal near the shadow edge has not previously 
been calculated by a time domain method. In this study, 
the response has been computed at all regions, from the 
illuminated zone to the deep shadow, using the Cagniard- 
de Hoop exact ray theory. An approximate model of plane 
homogeneous layers was obtained by applying an earth 


flattening transformation to the homogeneous earth model. 
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6.2 ‘The Homogeneous Earth Model 


The earth model shown in Figure 6.1 is not inten- 
ded to quantitatively satisfy experimental seismic data, 
but the model is simple and has been used in previous 
studies investigating other methods of seismogram cal- 
culation (Alexander and Phinney, 1966; Phinney and 
Cathles, 1969; Chapman and Phinney, 1972). Methods 
tested with this standard model can then be applied 
with confidence to realistic earth models. 

This model consists of a homogeneous liquid core 
in a homogeneous solid mantle. The dimensions corres- 
pond to the real earth with the crust removed. The 
parameters of the homogeneous mantle are those occur- 
ring at the core-mantle interface in the earth model 
used by Wiggins (1968), based on Birch (1964). All 
gradients have then been removed. In the resulting 
homogeneous earth model, all wave propagation paths 
are straight lines, and the geometrical shadow edge 
of the fluid core is located at i Si Wont Be, 
velocities and densities in Figure 6.1 are in km, 


km sec +, and gm cm? respectively. 


Orr The Pulse 


Smoothed delta function or inverse square root 


functions are unlike any features seen on real seismo- 
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HOMOGENEOUS EARTH 


MODEL 
mantle 
V,=-13.64 
V2 7.2 core 
rho=5.43 ea 
Le 


Figure 6.1 
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grams. Most seismologists prefer to convolve numerical 
seismograms with a wavelet function. To facilitate 
comparison, I have used, as a pulse, the step function 
response of a long period, critically-damped, coupled 
30 second seismometer and 100 second galvanometer as 
used by Chapman and Phinney (1972). The long period 
pulse is used because diffracted waves are low frequency 
Signals, detectable on long period instruments. The 
wave solution method is only applicable to band-limited 
pulses, and a low frequency pulse can be easily synthe- 
sized from a narrow frequency band. 

The pulse is defined by the transfer function 

2 


F(w) = ee a ee (6.2) 
(08 ieiek iwA) (we - Rew iws) 


Qo and W, are the natural frequencies of the galvano- 
meter and seismometer, and A = 220 and 6 = 205 Or 
critical damping. The shape of the time domain pulse 


is shown in Figure 6.2. 


6.4 The Cagniard-de Hoop Numerical Seismograms 


Since signal amplitude is the most important 
parameter of diffracted waves, kinematic or geometrical 
ray theory is inadequate for diffraction studies. The 


response solutions obtained by means of the Cagniard- 
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de Hoop exact ray theory, combined with a ray expansion 
and an earth flattening transformation, theoretically 
include the diffracted signal in a geometrical shadow. 
This has been confirmed in this study. Numerical, long 
period seismograms have been computed for the model in 
Section 6.2, ,from 76°, in the 71) luminated region, across 
the shadow edge at 113.8°, and into the deep shadow at 


13257 Or P, SV, and SH waves. 


6.4.1 Selection of the Layered Model and Ray Expansion 


After the earth flattening transformation of 
Figure 4.2 has been applied to the spherical model, the 
resulting flat model, with continuous gradients of elas- 
tic parameters, must be approximated by a stack of plane 
homogeneous layers. If the layers are not sufficiently 
thin, and individual rays can be resolved in the final 
seismogram, then the layered model cannot be considered 
equivalent to the original continuous model. The 
straightforward method to obtain a suitable model would 
bei tovstart with a very simple structure, consisting of 
a few thick layers, and calculate the total response by 
summing the response of all important rays. The order 
of a ray is the number of reflections it undergoes 
between the source and the receiver. The lower order 


rays are, in general, the most significant. It would 
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then be necessary to increase the number of layers in 
the model by several steps, at each step computing the 
complete response. When the complete response was no 
longer appreciably changed, the layered model would be 
equivalent to the continuous model. Because this 
process would require large computing times, it is 
sufficient to use the method described by Muller (1970). 
At each step, as the number of layers is increased, only 
the response of the first order rays is calculated. 
When this partial ray sum becomes stable, the partial 
sum due to the neglected higher order rays is also stable 
(Miller, 1970). With the sufficiently detailed model 
determined, higher order rays are computed and added to 
the ray sum until the total ray summation is stable. 
This process must be repeated over the range of epicen- 
tral distances of interest, because the nature of the 
required model and ray expansion is distance dependent. 

Layers corresponding to 40 kilometre shells in 
the spherical earth model were found to be sufficiently 
thin both near the turning depth of the direct wave in 
the illuminated zone, and near the core-mantle boundary 
at ranges corresponding to the core shadow of the 
spherical model. 

The time-sampling interval used for these seis- 
mograms was two seconds, and the response was smoothed 


by a Bartlett window (Figure 5.2) with a half-width of 
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four seconds. Shorter period seismograms would require 
thinner layers in the model to duplicate the response of 
the original continuous model. 

During this model selection process, it was found 
that the* partial ray sum for rays of a given order j 
(i.e. all rays in the partial sum have j reflection 
points) becomes only meta-stable with increasing number 
of layers. Figure 6.3 a) illustrates this using the SH 
diffracted signal deep in the shadow, at 124°. The 
Signals can be interpreted as either the time integral 
of displacement for the source potential in (3.20), or 
as the displacement response for a source displacement 
fUnCELOnN, DLOpPOEtiOnal to t7%, With layers corresponding 
to 80 kilometre thick spherical shells, the response due 
to first order rays has already stabilized. As the shell 
thickness is decreased to 10 kilometres, however, the 
response due to first order rays starts to diverge. 

The total ray summation appears to be truly 
stable, however. In Figure 6.3 b), third order rays 
Having two reflection points at the core have been 
added to the first order ray response. This response 
stays stable down to 20 kilometre shell thickness. The 
response from the 10 kilometre layered model is much 
Closer to that of the other models when these few third 
order rays are included, and the addition of more higher 


order rays will tend to reduce the difference further. 
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SH FIRST MOTION AT 124° 


IN HOMOGENEQUS EARTH MODEL 
TIME 

feeds on 1560. 
~~ THICKNESS OF LAYERS 


a) FIRST GRDER RAYS ONLY 


TIME 
Se 
bosu. 1560. 


b) FIRST ORDER RAYS AND SOME 
IMPGRTANT THIRD GRDER RAYS 


Figure 6.3 


Divergence of the First Order Ray Expansion for Thin Layers 


° ay Ret Aaa Ai 
i nt oy av ea : 
“usr TA MOLTeM TRA _ 


J300m HTAA 2UCaMIa0MOH ne 


: a : 
032 pu ORAL 
2ASYAS 48 BeSMNIIAT A ae 
MAX D8 i a 
ae Mx OH \ 
LA ad — wos 4 
& MNO! | eb 
‘ 
nb 


‘a 
- - 7 a. > 
: 7 | : 
fo : . 7 ; thks 
a - BAO? DHA SYAR ABOARD TAL @ 
7 a a Pa 7 eyAR ABNRO 1 isi bal HATO a 
7 Wie. fo 1 Ave ef 
| AYR ef 
7 - 7 x 


ae Cate Mi ee 
fy “, Sh. | S-elei2is .- vt 


80 


y°9 eanbta 


SONS tSSt eee dsa te SaavseNG 
INBIYOdWI Sddil Abe Ys0Y¥O YSHITH 


JAYS | | q ne 


TWRTAOSMT e977 TAR ASGRO ASHOTH 
eJAvdTe OSTSAAAITO-SABD UI 


bo stupr4 


81 


G*9 eanbtgq 


JABM LISYTO NI Sib Ys0u0 YSHOTH 
LNULYOdWI LSOW SHL INT LEY4NS9 


3 
Hi d30 
ONINENL 
Qq e 
Hida as H1d20 
SNINUAL ONINUNAL 


——— 


2 


MEAL? 


———— 
OMTMAGT — 


<a VA ag: 


TMATAOWML T28M SHT aALTAAIUSE 
SVAW TOIATG Wl eYAR ASOAD ABHOIH 


—§ > : . ¢.8 suupt’ | 7 


6&2 


When more layers than the sufficient number are 
used to model a continuous velocity gradient, it may 
be necessary to include proportionately more higher 
order rays in the ray expansion. Merely increasing 
the number of layers may actually decrease the accuracy 
of the computed response. This difficulty in ray ex- 
pansion studies is not widely appreciated. 

The relative importance of higher order rays 
might be expected to increase with the number of layers. 
The number of rays of order j in the complete ray 
expansion (3.17) increases like KI with increasing 
number of layers K. Higher order rays increase rapidly 
in abundance relative to first order (j=1) rays as K 
increases. 

When a model with a sufficient number of layers 
has been chosen, it is necessary to determine an appro- 
praate Layeexpansion.yeTosthel sum ofe first’ order rays 
are added third order rays, starting with those 
obviously most important. When the response is no 
longertsignificantly altered, the’ processyist repeated 
Withetwetheorder rays, fandgsoion: 

Thetmosite Lmportantecilrstsorderhraysinrthetdirect 
waveein the illuminated*region isrthationenreflected 
nearest the turning depth of the geometrical ray in 
the continuous model. The important higher order rays 


must carry energy close to the path of this ray. By 
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adding or subtracting pairs of ray segments, this ray 
can be modified to produce the most important higher 
order rays, as shown by Figure 6.5. In the illuminated 
region none of these higher order rays was significant. 

As the range was increased into the core shadow, 
the higher order rays became increasingly important. 

This behaviour of the ray expansion is easily 
understood. The path of energy in the direct wave in 
the illuminated region is very similar to the path of 
first order rays reflected in the layered model near 
the geometric turning depth. It is logical that these 
rays contain most of the energy in the direct wave. 

The energy in the diffracted signal in the core 
shadow has travelled a large distance near the core, 
and influenced by the core. No first order rays do 
this. Only one first order ray touches the core even 
once. Rays multiply-reflected in the lower mantle 
can travel large distances near the core, and are 
influenced by the core at a variety of reflection points 
along the core mantle boundary. These higher order rays 
are expected, on physical grounds, to contribute subs- 
tantially to the diffracted pulse. In a mathematical 
SOntext, it 1S evident, that at large distances; as ray 
segments near the core approach grazing incidence, the 


vertical wavenumbers n approach zero. The magnitude of 
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Sach reflection coelficient Six in Table 5.1 approaches 
unity. With the amplitude of a ray no longer depending 
strongly on the number of reflections, higher order rays 
become important. 

In the core shadow, four types of higher order 
tay, diavang all retlection points within the bottom 
ten layers of the mantle, were significant in the first 
23) Seconds, Of Ground motion. © Third order cays having 
two reflections at the core mantle interface (Figure 
6.4 a)) were the most important. Third order rays with 
one reflection at the core, and the other two reflection 
points separated by one layer (Figure 6.4 b)), or by two 
layers (Figure 6.4 c)) were also necessary. The diffrac- 
ted SV signal’in the deep shadow also required the 
similar third order rays reflected at the core 
and at two other interfaces separated by three layers. 
Other third order rays had arrival times or amplitudes 
outside the range of interest. 

(ie on yusclLonir cant clrtin order aay 1 tne mec 
shadow is shown in Figure 6.4 qd). 

No rays reflected within the core were important 
at the times and distances studied. The core was appro- 
ximated by a half space in the flattened model. This 
approximation would not be valid at larger ranges where 


core phases, such as PKP, appear, or at ranges where SKS 
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No rays with segments of converted wave type in 


the P and SV problems were significant. 


6.4.2 The SH Seismograms 


Long period numerical seismograms were computed 
Otecoelotervals, trom 76°. to. 132° for uthes model on 
Section 6.2. Figures 6.6 to 6.9 each show four stages 
in the SH ray summation. These figures cover ranges 
from the edge of the illuminated region, to deep in the 
core shadow. The signal is the displacement response 
to a source displacement of the form t7%, Line a) of 
each figure is the sum of first order rays only. Line 
b) includes the third order rays reflected twice at the 
core (examples in Figure 6.5 a)) and the fifth order ray 
Mier 1 OUT eCn6. 5.0)... bine ¢c) includes third-order rays Jot 
the type illustrated in Figure 6.5 b),, and line d) in- 
cludes third order rays of the type shown in Figure 6.5 
on 

The summation of these rays gives a good repre- 
sentation of the displacement for the first 20 seconds 
AG terotie arrival Ot, the Gifiracted P pulse, ana for 
30. seconds after the arrival.of the diffracted 5S pulses. 
The shape of the signal, as its amplitude decreases 
toward zero after 30 seconds, is undetermined, due to 
the use of a truncated ray expansion, and the Bessel 


function approximation (3.19). The final ray summation, 
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i.e. line d) in each figure, was truncated after 30 
seconds of motion and arbitrarily extended, at 
constant amplitude, to the end of the time window 
cae (Chapter 5). The complete response curve after 
the arrival of the diffracted pulse was then rotated 
to intersect the time axis 100 seconds after the 
diffracted pulse arrival. The response after 100 


seconds was then assumed to be zero. Figure 6.10 


shows these tapered long period SH seismograms for 


wv 


a source displacement of the form t *. The emplitudes 
are in metres, when multiplied by the source strength 
Fe in SI units. 

At 108° in the illuminated region, the signal 
has a shape similar to the source. Figure 6.6 shows 
that the third order rays are of very minor importance. 
As the range increases into the core shadow, the signal 
changes in character toward a broad, low frequency pulse. 
Figures 6.7 to 6.9 show that this is correlated with the 
increasing importance of the higher order rays. This 
change in pulse shape is expected. As in all diffrac- 
tion phenomena, the decay length of a diffracted wave 
increases with wavelength. For any given opaque body, 
long wavelength components form less distinct shadows 
than do shorter wavelengths. 

In order to compare these results with published 


seismograms obtained by a different method (Chapman and 
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Phinney, 1972), the signals an Figure 6.10 have been 


we 


deconvolved from the source displacement of t “, and 
convolved with the pulse of Section 6.3. This pulse 
was also used by Chapman and Phinney. The resulting 
Signals in Figure 6.11 may be interpreted as the 
displacement of the instrument: with transfer function 
(6.1), for an impulsive source displacement. Alter- 
natively, the signals are the displacement of a 
critically-damped, coupled 30 second seismometer and 
100 second galvanometer, for a Heaviside step source 
displacement. The rapid amplitude decay with distance 
into the shadow is best seen in this figure. 

Figure 6.12 shows the general behaviour of the 
response over the larger range of 76° to 132°. At 
76°, in the illuminated region, a direct wave, and 
the following smaller core reflection, are evident. 

As the shadow edge at 113.8° is approached, the time 
difference between these two arrivals goes to zero, and 
the decay rate for the single pulse in the shadow is 


much greater than for either signal in the illuminated 


region. 


Gets The P and SV Seismograms 


Similar computations for P waves are shown in 


Figures 6.13 to 6.19. The behaviour of the response 
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is very similar to the SH response discussed in the 
preceding section. 

Computations of long period SV seismograms are 
snown by “Figures 6.20 to 6526. The behaviour of the 
SV signals is fundamentally different from that of the 
PHeanouP signals. Ihe particle motion tor a Pon SH 
wave grazing the core is tangential: to the core boun- 
dary. “he Su owave 1s totally reflected by the fluid, 
and the grazing P wave also does not effectively transmit 
energy into the core. The particle motion ELor a grazing 
SV wave, however, is perpendicular to the core boundary. 
Une SV ysighal as altenuated much more rapidly than P or 
SH signals, as the range increases. 

It was necessary to include another set of higher 
SLdet Mays, OLier than those illustrated an figure 6.5, 
£o obtain the early diffracted SV motion an the score 
shadows Third older “ays with one reflection at the 
Core, ana two other reLlection points: separated by 
three layers. were found to alter the response. Line 
Sor Pi1guress6.201bO 6.25 cnowe tie -eLicCl One nciud— 
ing these rays. 

The SV signals in the shadow have a much dif- 
ferent shape than the P and SH signals, due to the 


existence of a P head wave along the core-mantle. 
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boundary. At ranges of 108° and beyond, the first 
arrival is this P head wave, and the later, shorter 
pulse, of the opposite polarity, is the diffracted SV 
Signal. 

The calculation of each ray in the SV ray 
expansicn is complicated by the existence of P head 
waves along each interface in the mantle of the 
layered flat model. The angle of the SV ray incident 
at an interface (see Figure 2.2 b)) is subcritical 
for SV waves below the interface, so a transmitted SV 
ray exists. However, the incident SV ray may be post- 
critical for P waves below the interface, creating a 
P head wave that is included in the exact ray response. 
This type of head wave cannot occur with incident P or 
SH waves. 

These P head waves do not correspond to any 
real wave in the original model, but have been intro- 
duced by the process of approximating that model by 
homogeneous layers. The computations showed that, at 
the ranges studied, these head waves are several orders 
of magnitude smaller than the main SV pulse, and tend 
to cancel one another. They do not significantly 
affect the final seismograms. They are significant 


mainly for their nuisance value. 
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Figure 6.22, showing the stages of the SV ray 
summation at 124°, includes, above line a), the res- 
ponse of only the single first order ray reflected 
from the core. This figure effectively illustrates 
the importance of the higher order rays in the core 
shadow. The first order ray above line a) contains most 
of the energy of the P head wave, but contains nothing 
at all similar to the final diffracted SV pulse in line 
6). ~The sumlot sfirst order rays) inline 4) "superiacially 
resembles the final line e), but the similarity is coin- 
cidental. The pulse resembling the diffracted signal 
disappears completely at line b) when the third order 
rays like Figure 6.5 a) are included. The diffracted 
pulse only appears as the other third order rays are 
added. Higher order rays are absolutely necessary in 
order to compute the diffracted signals. 

Thevaccuracy of the SV signals, at. 132° 1s quées—~ 
tionable. The ray expansion valid for 124° was used 
at 132°, and it is possible that even more third and 
fifth order rays are significant at 1322 - 

Early motions in the time domain response have 
been found accurately by the Cagniard-de Hoop ray 
method. There will, however, be uncertainties in the 
frequency spectrum of the total signal due to the 
necessity of arbitrarily tapering the response after 


30 seconds. 
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Ono Comparisons with the Wave Solutions 


Chapman and Phinney (1972) have published numeri- 
cal seismograms for P, SV, and SH signals over the same 
epicentral distances in the same earth model used in 
this study. Using the wave method outlined in Chapter 1, 
Chapman calculated the time-transformed response at 
sixteen frequencies from .0156 Hz to 0.25 Hz. Sixteen- 
EOtdyirequency interpolation, and multiplicaticnrby 
the pulse transfer function (6.1), gave 256 spectral 
values. An inverse Fourier transform produced the 
seismograms of the displacement potential for SH, P, 
and SV waves ,; as: shown jin Figures 6227 ito 6.29 (irom 
Chapman and Phinney, 1972). The normalisation used is 
not important here. 

These potentials show the same qualitative be- 
haviour as the Cagniard-de Hoop displacements of 
Figures 6.121 6.19, wand 6.26. ° Aldirect pulse sand a 
weaker, later’ reflection from’ the core are resolved 
in the illuminated region. These two pulses merge 
Near the shadow edge ati b1ls.8°, andthe single pulse 
is highly attenuated with distance into the shadow, 

SV more So than P Or sk. 

A quantitative comparison of pulse shapes and 

of signal decay rates is possible in the frequency 


domain. Each solid curve in Figures 6.30 to 6.35 shows 
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the distance dependence of the amplitude of the complex 
spectrum of the displacement response at one frequency, 
as calculated by Chapman. The dimensionless frequency 


parameter ka is defined by 


ka = wa/a (622) 


a® is the P wave velocity in the mantle, and a is the 
outer radius of the earth (see Figure 6.1). The 


i eee aueneced 


curves have been normalized to (kR) 
region. R is the path length of the direct ray. The 
crosses and circles are the corresponding spectral 
amplitudes of the Cagniard-de Hoop seismograms, also 
normalized to Cen) 

The wave spectral amplitudes in the illuminated 
region show beating due to the two close pulses in the 
time domain. The decay of each spectral component is 
approximately exponential across the shadow. The 
decay constant, given by the inverse slope in these 
logarithmic plots, decreases with increasing frequency. 
This results in a relatively low frequency pulse in the 
deep shadow as was noted in Section 6.4. 

The correspondence of the wave and the ray solu- 
tions is quite acceptable when the differences in com- 
putational methods are considered. At the lower fre- 


quencies, the amplitude values consistently agree within 


30%. At higher frequencies, the agreement is not as 
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SH Spectral Amplitude Comparisons, Lines are Wave Solutions. 


Crosses and Circles are Ray Solutions. 
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good. Aliasing in the Fast Fourier transformation of 
the ray solutions is not important below ka = 200, but 
the nature of a ray expansion, as a summation of spikes, 
Pamits ie accuracy at high frequencies. 

The amplitudes in the illuminated region are 
generally larger for the ray seismograms than for the 
wave solutions. The linear interpolation of the con- 


-% 
cave t 


response, while smoothing and digitizing each 
ray (Chapter 5), results in an over-estimation of the 
amplitude. The response can be improved by a quadratic 
interpolation method. 

The decay rate of the Cagniard-de Hoop spectral 
amplitudes tends to be greater than that of the wave 
solutions. The decay constants in Table 6.1 are the 
angular range, in radians, over which the spectral 
amplitude decreases by l/fe. These values have been 
determined by fitting the best straight line to the 
wave solution and ray solution spectral amplitudes 
across the shadow boundary from 108° to 124° on the 
logarithmic plots an Figures 6.30 °t6 6.35. “Thersv ray 
solution data at ka = 150 and 200, were not close enough 
tova ‘straight Dine*to be useful. © Because the error’ in 
the ray solution decay constants tends to decrease with 
increasing frequency ka, the discrepancy with the wave 


results appears to be due to the arbitrary method of 
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tapering the ray seismograms after 30 seconds of 
motion. The response pulse lengthens with distance 
into the shadow, and tapering the response to zero 
after 100 seconds has filtered out the low frequency 
energy that should arrive later, in the broad diffrac- 
ted pulse. The ray solutions were designed to give 

a total response at all frequencies in a time window 
of the first 30 seconds of motion, while each wave 
solutior is specifically designed to find an accurate 
amplitude of a single frequency component in the com- 
plete pulse. A large part of the discrepancy in decay 
rates predicted by the two solutions is introduced 
Ehrough the basic, conceptual difference between wave 
and ray methods. 

The poor agreement of spectral amplitudes for 
the OV signals fat 32° confirma that more naghes order 
rays are necessary in that seismogram. 

Considering these limitations, it is concluded 
that acceptable agreement can be obtained between wave 
solutions and Cagniard-de Hoop ray solutions, under 
conditicns extremely favourable to the wave method. 

Models containing variable velocity gradients, 
or thin, resonating layers, have more complex frequency 
spectra. Interpolation between wave solution spectral 


amplitudes is very difficult. Ray solutions, however, 
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are no more difficult for such a radially inhomogeneous 
model™than* for ’thevsimple model -of- this study. 

The Cagniard-de Hoop exact ray method appears to 
be a practical way to compute numerical seismograms for 
studies of the core-mantle boundary, and will be parti- 
cularly useful for models with rapid radial variations 


of elastic parameters. 
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CHAPTER 7 


NUMERICAL SEISMOGRAMS FOR A STUDY OF LOWER MANTLE 


HETEROGENEITY 


Jv leySetremic Arrays 


The introduction in recent years of coordinated 
arrays of seismometers, spread out over distances of 
tens to hundreds of kilometres, has enabled seismologists 
tO ‘accurately measure the azimuth of arrival and the 
phase velocity dA/dT of seismic waves. A is the dis- 
tance from the earthquake epicentre to the array, and 
T is time. This information is a valuable supplement 
to data on group velocities, determined by arrival times 
at individual stations. The inverse phase velocity, or 
wave slowness dT/dA, is used in the Wiechert-Herglotz 
mniversrontprocess ((Bulleny 1963; op 9119) seeeobtainathe 
velocity-depth function of a spherically symmetric earth 
from travel time data. The slope of the travel time 
eunvesiars thertraditionalasource vof -dT/dA *data, abut 
direct measurement of dT/dA by seismic arrays is reduc- 
ing uncertainties in the travel time curves and in earth 
Struceucer 

The Large Aperture Seismic Array (LASA) in 
Montana, having dimensions of the order of 300 kilometres, 
was built to distinguish nuclear weapons tests from 


natural earthquakes. 
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The smaller Yellowknife array in Canada, and 
the Warramunga array in Australia have dimensions of 
the order of 20 kilometres. 

The Variable Aperture Seismic Array (VASA) is an 
array operated by the University of Alberta exclusively 
for purposes of geophysical research. It consists of nine 
broad band (.03 to 4.0 Hz), three component seismographs 


mounted in mobile trailers, in arrays of intermediate size. 


7.2 Lower Mantle Inhomogeneity Beneath Hawail 


The VASA, located in central Alberta, has measured 
aT/dA for earthquakes at a wide distribution of azimuths. 
Most values agree well with the slopes of the Jeffreys- 
Bullen travel time curves (Kanasewich et al, 1972). 
However, waves from earthquakes at Tonga and Samoa, 
at Bercener at distances of 84° to 95°, have extraordinarily 
high phase velocities. Figure 7.1, from Kanasewich et al, 
1972, shows that events in the Solomon Islands, Asia, 
and South America have phase velocities (shown as crosses) 
close to the phase velocity (the solid line) expected in 
the Barth Model IL of Burehe(1964). The phase velocities 
of the Tonga and Samoa events are shown as triangles. 

Kanasewich et al (1972) have shown that the anoma- 
lies are probably due to inhomogeneities at the bottoming 


point on the energy propagation path. Geometric ray 
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MODIFICATION OF 
BIRCH EARTH MODEL O 
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Hawaii Model (Solid Lines) and Birch Earth Model II 


(Dashed Line). 
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theory indicates that energy arriving at VASA from Tonga 
and Samoa bottoms within several hundred kilometres of 
the core, beneath an area several hundred kilometres 
south-east of Hawaii. Kanasewich et al (1973) 
have proposed modifications of the Birch Earth Model II 
in the vicinity of Hawaii. This Hawaii model is shown 
by the solid lines in Figure 7.2. The original model 
due to Birch is illustrated by the dashed line. 
Kanasewich et al suggest that the phase velocities 
in Figure 7.1 that are 2 km/sec higher than expected 
between 84° and 90°, are due to an increased P wave 
velocity gradient between 4000 and 3500 kilometres radius. 
The sharp jump to phase velocities over 3 km/sec greater 
than expected may be due to propagation in a layer 
approximately 30 kilometres thick at the base of the 
mantle. This waveguide is postulated to have a P wave 
velocity approximately 0.7 km/sec greater than that in 
the mantle above. 

Observations of low phase velocities and anomalies 
in azimuth for some events recorded at LASA (Davies and 
Sheppard, 1972) are helping to delineate the lateral 


extent of the inhomogeneities beneath Hawaii. 
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AS NS Cagniarcd-de Hoop Seismograms 


A successful model for earth structure must 
predict correct amplitudes and pulse shape, as well 
as the purely kinematic phase velocity. Numerical 
seismograms for models having thin layers, with large 
velocity contrasts, are difficult with wave methods. 
Interference of waves in the thin layers causes reso- 
nance in the frequency spectrum. The interpolation of 
intermediate frequencies cannot be done with any degree 
Of certainty; SO it is necessary to explifeitiy calculate | 
Many more spectral amplitudes than were required for the 
homogeneous model in Chapter 6. 

Ray methods are not limited in this way. A few 
Pays, multipiy-retlected ian a high velocity Layer, can 
be added to the ray summation quite easily. The 
Cagniard-de Hoop ray theory is the method best-suited 
to the study of a thin layer problem near the geometric 
shadow boundary. 

Numerical seismograms for epicentral distances of 
84° to 96° have been calculated for the Birch Earth 
Model II, and for a suite of models all having the 
increased P wave velocity gradient shown in Figure 7.2, 
but differing in the parameters of the thin layer at the 
base of the mantle. These models are defined in Figure 


7.3. Models I, II and III have high velocity waveguides 
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of varying thickness, bracketing the value of 27 
kilometres suggested by Kanasewich et al (1973). 

Models IV and V_ are designed to reveal the 

effects of changing the velocity contrast across the 

top interface of the thin layer. Model IV has a velo- 
city reversal, as has been suggested in a variety of 
contexts by several authors (Phinney and Alexander, 

1966 and 1969; Davies and Sheppard, 1972). Model V 
illustrates the effect of the increased velocity gradient 
alone. 

The pulse used in this study was the step func- 
tion response of a short period, critically-damped, 
coupled 1 second seismometer and 15 second galvanometer 
palr. The tyansfer function has the form (6.1), wath 
Ue and ips as 1 and 15 seconds, and A = 200 and 6 = 206° 
The time domain pulse is shown in Figure 7.4. 

The earth flattening transformation of Figure 
4.2 was applied to each model, and appropriate layer 
thicknesses and ray expansions chosen as in Chapter 6. 
Because a short period response is required, it was 
necessary to use thin layers, corresponding to shells 
25 kilometres thick near the core-mantle boundary in the 
spherical model. At all ranges, only first order rays 
were necessary in all sections of the model above the 


thin waveguide layer. For the modified models I, II, IIil 


and IV, third order and fifth order rays, multiply- 
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reflected inside the thin layer, have been added to the 
ray sum. Rays with one segment traversing the wave- 
guide as a shear wave have been included. 

The }displacement, at ranges of 84° to 96° in 


the normal Birch Model II, for a displacement source 


nr 


proportional to t *, is shown in Figure 7.5. Displace- 
ment units are metres, when multiplied by the source 
strength Fo Hysol sun ts, . has wisethe basic. form of sine 
seismograms computed by the Cagniard-de Hoop method. 


-L 
* source has been deconvolved, and 


im Figure j26, the t 
the result convolved with the pulse in Figure 7.4. 

Since this normal earth model has a geometric core 

shadow in the vicinity of 96°, the seismograms in 

Figure 7.6 are in the illuminated zone, up to the shadow 
boundary. The direct P wave, and the core reflection 
(marked by the arrow), are clearly resolved at 84°. 

The two pulses merge by 90° and the amplitude decreases 
sharply from 93° to 96°, as the shadow edge is approached. 
This behaviour is expected by analogy with the diffrac- 
tion studies in Chapter 6. 

The minor ripples in the seismograms after the 
reflected signal would not be present if more later- 
arriving rays were included in the ray expansion. 

The displacement response of the modified models 
I through V is shown by Figures 7.7 to 7.11. Examina- 


tion of Figure 7.11 for model V shows that the inclusion 
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of only the increased P wave velocity gradient has 
moved the geometrical shadow edge closer to the source, 
so that the range 84° to 96° spans the shadow edge. 
The direct and reflected pulses are not resolved, and 
the amplitude decays steeply with increasing range. 

Pi'gures*/.#acoeds9 show that adding a high velo- 
City waveguide at the base of the mantle creates much 
more structure in the seismic signals. The most pro- 
Minent new feature is the low frequency head wave along 
the top of the high velocity layer. Due to the higher 
velocity of propagation along this layer, the head wave 
has a higher dA/dT than, and arrives progressively earlier 
relative to, the main P wave. This is also true of the 
P wave (POP), which penetrates the high velocity layer 
and—rvs-reftbected by the core. Beyond 87°, the first 
motion is this short, weak POP pulse (marked by arrows), 
followed immediately by the broad, emerging head wave. 
The much larger P phase arrives later. The main P phase 
issalways the’ first arrival in the normal Birch model I1 
(Higuye=/..6),7 andiin ithe modets fiV and V (Figures 7.10 
and 7.ll) without @ high velocity layer. A weak precur= 
sor wave is apparently diagnostic of a high velocity 
layer near the turning depth. 

Figure 7.12 a) shows’ the record of an event from 
the Solomon Islands that had a normal phase velocity. 
The first arrival is the P wave, as in the normal earth 


model (Figure 7.6). Figure 7.12 b) is an event from 
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a NORMAL EVENT 
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b ANOMALOUS EVENT 
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Tonga, having a high phase velocity. The weak precursor 
wave could be a head wave on a high velocity layer, and 
the PUP phase. 

A direct comparison of pulse shape and Spi ede 
is not feasible at this stage, because the earthquake 
source functions are unknown, and cannot be deconvolved 
from, thesdeta records in) Figure 7212. The numerical 
results assume a step function source displacement. 

The receiver functions are also different, although this 
could be corrected by the appropriate convolution opera- 
tions. It is appropriate to note, however, that the 
relative amplitudes of the precursor and the P wave of 
the numerical seismograms of Figure 7.8 for the modified 
model II are compatible with the observed amplitudes. 

For model I, having the thinnest high velocity 
layer, the second arrow in Figure 7.7 indicates energy 
that has travelled once through the layer as a P wave 
and once as an SV wave, and as a P wave both ways in 
the mantle. The importance of these rays with converted 
segments is expected to increase as the ray approaches 
grazing incidence at the core, due either to increasing 
range in a-“fixed layer, or to décreasing Payer thickness 
ata fixed range. When the vertical wavenumber n ap- 


proaches: zero), theodifference’an° magnitude of the  P-P 
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and P-SV reflection coefficients § and S (Table 


le 2 

5.1) disappears, as they both approach unity. This 
weak pulse would possibly be visible on seismic records, 
under good signal-to-noise conditions, if a very thin 
high velocity layer is present. 

The response in model Iv, with a low velocity 
layer, is quite different. Three phases are evident. 
The P wave is always the first arrival. This is 
followed by a smaller pulse of the opposite sign 
(first arrow in Figure 7.10), reflected at the top 
of the low velocity layer. The third arrival, marked 
by the second arrow, is the core reflection. Due to 
propagation in the low velocity layer, this pulse 
arrives significantly later than the P wave. “The 
existence of two separate pulses, even in the P wave 
shadow, as defined by the exponential amplitude decay 


rate, appears to be diagnostic of a low velocity zone 


at the base of the mantle. 


72s Conclusrons 


The question of the validity of modelling a 
region of known lateral heterogeneity by means of 
homogeneous layers of unbounded extent, is an impor- 


tant one. The procedure is expected to be valid if 
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the propagation path length in the region is much 
Shorter than the scale length of the inhomogeneities. 
The dimensions of the high velocity anomaly under 
Hawaii appear to be several hundred kilometres - 
(Gutowski, personal communication),which is less than, 
or comparable to the propagation path length below 
4000 kilometres radius. Gutowski has attempted a 
more quantitative comparison of these numerical seis- 
mograms with the observational data, by means of a 
spectral ratio technique, designed to remove source 
effects. His preliminary results (personal communica- 
tion) indicate that the structure at the core-mantle 
boundary in Birch's Model II (1964) may be inadequate 
for the numerical duplication of normal seismic events 
such as in Figure 6.12 a). 

However, with these two limitations in mind, 
the qualitative agreement of the numerical seismograms 
of the modified model II (Figure 7.8’ with the observed 
data» (Figures “7.1 and 7-22 5)) as encouraging. The 
main features of the model proposed by Kanasewich 
Sb eas (bo 3) appear to be correct. Although wave 
methods are inappropriate for constructing numerical 
seismograms for structural problems with thin resonat- 
ing layers, Cagniard-de Hoop seismograms will. be a use- 


ful complement to the kinematic array observations, in 
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the study of structure near the core mantle interface, 
especially for ranges near the core shadow boundary. 
Amplitude data from stations over distances of 10° can 
be used to recognize shifts in the shadow nye eee due 
to. altered, velocity gradients,.and to confirm, and 
define high and low velocity zones. 

The integrated use of the Cagniard-de Hoop exact 
ray method with the asymptotic ray method may make it 
feasible to include lateral variations in the earth 
models and help resolve many of these important earth 


structural problems. 
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CONCLUSIONS 


The response of earth models to elastic waves 
at the vicinity of the geometrical shadow of the liquid 
core has been studied by means of an earth flattening 
transformation, a ray expansion, and the Cagniard-de 
Hoop exact ray theory. The solution of this problem 
has not previously been obtained by ray methods. 

A sum of first order rays adequately represents 
the total signal in the illuminated Zone, but third 
and higher order rays are an essential part of the ray 
summation in the core shadow. 

The validity of the method has been verified by 
obtaining diffracted signals comparable to those pu- 
blished by Chapman and Phinney (1972) for a homogeneous 
core in a homogeneous mantle. Acceptable agreement was 
found between the early motion obtained by the exact 
ray theory with the approximate layered model, and the 
frequency band-limited signal obtained by an approximate 
numerical wave method with the exact model. 

The Cagniard-de Hoop ray method was then used to 
compute numerical seismograms for a proposed earth model 
having a structure unsuitable for wave method solutions, 
and at ranges unsuitable for previous ray methods. 
Qualitative agreement with observational data was 


obtained, and the results indicate that two and three 
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dimensional modelling methods must be developed to 
study recently discovered lateral variations in the 
lower mantle by means of numerical seismograms. This 
may be achieved by a combination of the Cagniard-de 
Hoop exact ray method to evaluate grazing and critically- 
reflected rays, and the asymptotic ray method to evaluate 
rays at curvilinear interfaces and lateral velocity 
gradients. 

The normal structure found at the core-mantle 
interface is still not well determined. The Cagniard- 
de Hoop exact ray method is most appropriate for inves- 
tigating. thin layers or, velocity reversals in the, dower 


mantle by comparing data with numerical seismograms. 
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